pacman::p_load(tmap, sf, sp,
performance, reshape2,
ggpubr, tidyverse,
DT, stplanr)Hands-on Exercise 11: Calibrating Spatial Interaction Models with R
1.0 Introduction
1.1 Getting Started
In this hands-on exercise, we will learn how to calibrate SIM to determine factors affecting the public bus passenger flows during the morning peak in Singapore.
1.2 Overview
Spatial Interaction Models (SIMs) are mathematical models for estimating flows between spatial entities developed by Alan Wilson in the late 1960s and early 1970, with considerable uptake and refinement for transport modelling.
There are four main types of traditional SIMs:
Unconstrained
Production-constrained
Attraction-constrained
Doubly-constrained
1.3 Installing and loading R packages
2.0 Data Acquisition
We will be using 2 datasets in this exercise:
Weekday morning peak passenger flows at planning subzone level (i.e.
od_data.rds)URA Master Plan 2019 Planning Subzone boundary in simple feature tibble data frame format (i.e.
mpsz.rds)
3.0 Geospatial and Aspatial Data Handling
We will be using the st_read() from sf package to import the data into RStudio.
3.1 Importing Geospatial Data
mpsz <- st_read("data/geospatial",
layer = "MPSZ-2019")Reading layer `MPSZ-2019' from data source
`C:\kt526\IS415-GAA\Hands-on_Ex\Hands-on_Ex11\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS: WGS 84
mpsz_sp <- as(mpsz, "Spatial")
mpsz_spclass : SpatialPolygonsDataFrame
features : 332
extent : 103.6057, 104.0885, 1.158699, 1.470775 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
variables : 6
names : SUBZONE_N, SUBZONE_C, PLN_AREA_N, PLN_AREA_C, REGION_N, REGION_C
min values : ADMIRALTY, AMSZ01, ANG MO KIO, AM, CENTRAL REGION, CR
max values : YUNNAN, YSSZ09, YISHUN, YS, WEST REGION, WR
Computing the distance matrix
dist <- spDists(mpsz_sp,
longlat = FALSE)
sz_names <- mpsz$SUBZONE_C
colnames(dist) <- paste0(sz_names)
rownames(dist) <- paste0(sz_names)3.2 Importing Aspatial Data
population <- read.csv("data/aspatial/pop.csv")4.0 Preparing Flow Data
4.1 Importing the Origin Destination data
odbus <- read_csv("data/aspatial/origin_destination_bus_202210.csv")
odbus$ORIGIN_PT_CODE <- as.factor(odbus$ORIGIN_PT_CODE)
odbus$DESTINATION_PT_CODE <- as.factor(odbus$DESTINATION_PT_CODE) 4.1.1 Extracting the study data
odbus6_9 <- odbus %>%
filter(DAY_TYPE == "WEEKDAY") %>%
filter(TIME_PER_HOUR >= 6 &
TIME_PER_HOUR <= 9) %>%
group_by(ORIGIN_PT_CODE,
DESTINATION_PT_CODE) %>%
summarise(TRIPS = sum(TOTAL_TRIPS))
datatable(odbus6_9)Writing to rds
write_rds(odbus6_9, "data/rds/odbus6_9.rds")4.2 Importing bus stop data
busstop <- st_read(dsn = "data/geospatial",
layer = "BusStop") %>%
st_transform(crs = 3414)Reading layer `BusStop' from data source
`C:\kt526\IS415-GAA\Hands-on_Ex\Hands-on_Ex11\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 5159 features and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 3970.122 ymin: 26482.1 xmax: 48284.56 ymax: 52983.82
Projected CRS: SVY21
4.3 Importing MPSZ-2019
mpsz <- st_read(dsn = "data/geospatial",
layer = "MPSZ-2019") %>%
st_transform(crs = 3414)Reading layer `MPSZ-2019' from data source
`C:\kt526\IS415-GAA\Hands-on_Ex\Hands-on_Ex11\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS: WGS 84
Write data to rds
mpsz <- write_rds(mpsz, "data/rds/mpsz.rds")4.4 Combing bus stop and mpsz
busstop_mpsz <- st_intersection(busstop, mpsz) %>%
select(BUS_STOP_N, SUBZONE_C) %>%
st_drop_geometry()
datatable(busstop_mpsz)Write data to rds
write_rds(busstop_mpsz, "data/rds/busstop_mpsz.rds") od_data <- left_join(odbus6_9 , busstop_mpsz,
by = c("ORIGIN_PT_CODE" = "BUS_STOP_N")) %>%
rename(ORIGIN_BS = ORIGIN_PT_CODE,
ORIGIN_SZ = SUBZONE_C,
DESTIN_BS = DESTINATION_PT_CODE)To check for duplicates
duplicate <- od_data %>%
group_by_all() %>%
filter(n()>1) %>%
ungroup()
od_data <- unique(od_data)
od_data <- left_join(od_data , busstop_mpsz,
by = c("DESTIN_BS" = "BUS_STOP_N")) duplicate <- od_data %>%
group_by_all() %>%
filter(n()>1) %>%
ungroup()
od_data <- unique(od_data)